Technical Report on the NASA Research Grant NAGS- 13388 "Validation of mean absolute 
sea level of the North Atlantic obtained from drifter, altimetry and wind data" conducted at 
the BPRC/SOEST, University of Hawaii, 05/01/03-09/30/03 


1. Introduction to the method 

Mean absolute sea level reflects the deviation of the ocean surface from geoid due to the ocean 
currents and is an important characteristic of the dynamical state of the ocean. Values of its 
spatial variations (order of 1 m) are generally much smaller than deviations of the geoid shape 
from ellipsoid (order of 100 m) that makes the derivation of the absolute mean sea level a 
difficult task for gravity and satellite altimetry observations. 

Technique used by Niiler et al. (2003a) for computation of the absolute mean sea level in the 
Kuroshio Extension was then developed into more general method and applied by Niiler et al. 
(2003b) to the global ocean (Fig.l). The method is based on the consideration of balance of 
horizontal momentum in its simplest form: 

dV/dt + f X V = - g V h -b 3T / dz /po + H.O.T., (1) 

where f is the Coriolis parameter and g is a gravitational constant. After neglecting the higher 
order terms (H.O.T.), horizontal gradient of the mean sea level can be estimated as 

V<h> = -(dV/dt + f X V )/g - V h’ + 0T / dz /(pog) (2) 

I n m IV 

Acceleration (I) and Coriolis (II) terms are estimated using velocities derived from Lagrangian 
drifter data received from the NOAA/AOML, where they were acquired, quality controlled and 
interpolated onto 6-hourly intervals. Subtraction of the anomalous sea level gradient (IE) 
corrects the bias in the drifter ensemble due to its nonuniform temporal distribution and 
interannual variations of surface circulation. Term (IV) describing vertical divergence of the 
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Ekman stress was estimated using NCAR/NCEP reanalysis wind data at 10 m level and 
parameterization formula suggested by Ralph and Niiler (1999; thereafter RN99). Originally 
derived for Ekman velocities, it can be rewritten for (TV) in the complex notation as 


aT/az = Of"^Wexp{i(Po}, 


( 3 ) 


where W is the wind vector and the angle (po = 36® sign(6) is constant for each of the two 
hemispheres. 

Among advantages of the method are clear definition of the "mean" as a time average over 
specific time periods and high, mesoscale, spatial resolution. Such resolution helps to reveal the 
complexity of the mean surface circulation in many regions. Example of the North Atlantic is 
shown in Figure 2. 

Largest errors in (1) are expected in the term IV. These errors are both due to the uncertainties in 
the wind data and parameterization formula. Additional errors occur in the sea level anomaly 
(HI) field, generally, oversmoothed by the Aviso mapping procedure. 

This document reports the results of the study of these three sources of uncertainties. For the 
convenience of the study the geographic area was expanded to include all regions covered with 
the data. Using recently released preliminary data of the GRACE project, we were able to 
improve the Ekman parameterization suggested by Ralf and Niiler (1999). Improvement is 
particularly significant at high latitudes. By comparing between the Ekman parameterizations 
obtained for the NCAR/NCEP reanalysis winds and for the QuikSCAT satellite winds, we show 
that results are sensitive to the product and suggest future multi-variable analysis. Section 4 
describes biases that may occur in the data statistics when too simple methods are used to 
eliminate the difference in spatiotemporal spectra of drifter, satellite altimetry and wind data 
caused by different resolutions of the three datasets. 
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2. Ekman parameterization 

Among many sources of the potential errors in the mean sea level computed by integrating (2), 
the most important one is expected to be related to the way how Ekman currents (or stress) are 
parameterized. RN99 derived their formula (3) for low latitude regions of relatively low eddy 
activity and steady winds in the North and South Pacific. Their formula gives optimal latitude- 
dependent angle and coefficient between the mean NCAR/NCEP wind and Ekman currents, the 
latter being estimated using mean drifter velocities and climatological mean geostrophic 
velocities referenced to 1000 m depth. Figure 3 shows that in terms of horizontal curl of vertical 
divergence of the Ekman stress the formula gives consistent estimates in all ocean basins at least 
between 40°S and 40°N. 

However, it is not clear if: 

- the same formula can be applied to the high latitudes; 

- coefficients are not very sensitive to the frequency of the signal; 

- relation between the wind and Ekman velocity/stress is linear at all. 

Questions on the accuracy of (3) also remain because of vague definitions of the mean (“drifter- 
ensemble-mean” and “Levitus-climatology-mean”) used by RN99. Being weaker than typical 
ocean currents, Ekman currents, however, can give rise to large errors in the sea level when 
integrated over large distances. 

Recently, first results from the GRACE satellite mission provided improved model of the Earth’s 
gravity field (http://www.csr.utexas.edu/grace/). This geoid model was then used to compute the 
mean absolute sea level from accumulated during the last decade satellite altimetry. The products 
developed in various institutions are summarized and intercompared at 
http://iprc.soest.hawaii.edu/~nikolai/Globe/GRACE/GRACE.html. For this study we used 
products from the NASA Goddard Space Flight Center (courtesy of Brian Beckley, Figure 4) 
and Jet Propulsion Laboratory (courtesy of Victor Zlotnicky, Figure 5). Figures 4 and 5 illustrate 
the fact that effective spatial resolution of the GRACE data is limited to 400-500 km. However, 
larger scales seem to be trustworthy. Figure 6 shows that biggest difference between the GRACE 
(Fig.5) and NMM03 (Fig.l) mean sea levels is in the Southern ocean and reaches 130 cm across 
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the Antarctic Circumpolar Current (ACC). It was also found that zonal gradients of the mean sea 
level estimated from (2) with RN99 do not vanish while integrated along quasi-zonally contour 
around Antarctida. Unbalanced sea level gradient along the ACC also results in average 
differences between the oceans seen in Figure 6. Other differences can be attributed to different 
scales of spatial resolution (most obvious around the Gulf Stream, Kuroshio and Agulhas jets) 
and effects of coastal friction and upwelling. 

As a first step we use (2) to estimate T = 3T / 3z /(Po) at drifter locations and times. Term 

V<h> is estimated from the GRACE-based mean sea level (Fig.5), terms (I) and (II) come from 
drifters and term (HI) is from the Aviso/Enact grided sea level anomaly data. Technically, the 
procedure was designed to eliminate high-level noise in the data without oversmoothing the 
results. It is easy to show that coefficients of a simple rotate-and-rescale transform 
Tx = a-Wx-bWy, 

Ty = b-Wx + a-Wy, (4) 

minimize «<(T- A W)^»> (W is the wind vector) when 
a=«<TW»>/«<WW»>, 

b=<«kx(TxW)>»/«<W-W>», (5) 

where k is an upward looking unity vector. Then optimal angle and coefficient between vectors 
T and W are: 

tan(angle) = b/a, 

coefficient = sqrt(aVb^). (6) 

Thus defined coefficients take into account both correspondence between mean vectors and their 
covariance. 

Averaging <...> has been done in a number of subsequent steps. First, all instantaneous data 
were sorted according to local latitude and instantaneous wind speed and averaged (<...>) within 
the bins 1® latitude x 1 m/s. Second, values of <T-W>, <kx(TxW)>, <W W>, <T’-W’>, 
<W’ W’> and <T’ T’>, were additionally smoothed with the running 9® latitude x 3 m/s mean 
(«...») and correlation coefficient was calculated as CC=«T’'W’»/(«W’-W’» • 
«T’ T’»)*^^. Finally, at each latitude «T W», «kx(TxW)» and «W-W» were 
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averaged over the bins containing more than 100 drifter fixes (that defined extreme latitude 
limits) and having CC > 0.1. These values were than substituted into (5) and (6) and optimal 
angle and coefficient of T to W are shown in Figure 7. 

It is seen that qualitative correspondence with RN99 is fairly good. For some unclear reasons 
correspondence in the angle is especially good at 20-50*^ and correspondence in the coefficient 
is especially good at 10-50°S. Qualitative correspondence with RN99 seems to be even better in 
terms of Ekman velocities (Fig. 8). Indeed, jump in the angle across the equator looks more 
significant than its variations away from the equator. Coefficient differs from RN99 by less than 
30% between 60°S and 70°N (except for the narrow equatorial belt where the use of geostrophy 
is tricky). However, after more careful look, in both hemispheres Ekman velocities are almost 
perpendicular to the wind at low latitudes but become more in the direction of the wind around 
60° latitude. Such a monotonic behavior of the angle can be explained by the increase of the 
Ekman layer thickness with the growing latitude. The thicker the Ekman layer the closer the 
Ekman currents measured by drifters at 15 m depth to their values at the sea surface. 
Interestingly, the 30° angle suggested at 60° latitudes cannot be explained by the Ekman layer 
model with any constant mixing coefficient, but it does can exist in the model where mixing 
coefficient grows with the depth. Light blue lines in Figures 7 and 8 are drawn by reversing the 
latitude axis. In terms of the angle correspondence between the two hemispheres is very good 
(two blue lines deviate not more than few degrees from each other). One can see that even in 
terms of T (Fig.7) the angle changes sharply by about 30° across the equator. We assume that 
equatorial dynamics can involve higher order terms neglected here. Simple rule is that the angle 
between T and wind is approximately equal to the value of local latitude. 

The following tendencies are seen in the coefficients in Figs. 7 and 8. New coefficients are 
smaller than suggested by RN99 at 30°S and 30°N, latitudes of the weakest winds. New 
coefficients are larger than suggested by RN99 at the equator where winds are strong and steady 
and at 60°S and 60°N, where atmospheric eddies are very active, especially, the winter storms. 
One explanation can be that actual relation is nonlinear to the wind speed. Among various 
possible sources of non-linearity most credible are the actual relation between the wind vector 
and the wind stress and influence of the wind strength on the intensity of mixing in the upper 
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ocean. When the former fact is well known from the bulk formula the latter one is illustrated by 
Figure 9, which shows that in the northem/southem hemisphere T (and also Ve) at 15 m rotates 
counterclockwise/clockwise as wind becomes higher. Such a behavior is consistent with the 
suggestion that higher wind enhances the mixing, which increases the effective Ekman scale and 
changes values of the Ekman spiral at 15 m toward those at the sea surface. This process can also 
be responsible for the larger coefficients (lower panel in Fig.9) at higher winds. More details on 
the dependence of T on the wind strength are shown in Figure 10 and discussed at 
http://iprc.soest.hawaii.edu/~nikolai/Globe/GRACE/OptEkman-GRACE.html. 

Remarkable in the meridional distribution of coefficients are maxima at 65®S and 55®N. More 
extensive study is required to understand what makes these high latitudes so special and why the 
peak in the south is much more pronounced than in the north. Possible mechanisms can be 
related to the air-sea heat fluxes contributing to the mixing and also the effect of the Stokes’ drift 
in the wind waves field (McWilliams and Restrepo, 1999). 


3. Sensitivity of Ekman parameterization to the wind product: QuikSCAT vs NCEP 

Procedure similar to the one described in the previous section was applied to the data of 
QuikSCAT satellite wind data. The data (Level 2 sea wind vectors) are on spatial 25 km grid and 
were downloaded by Dr. Jan Hafner (IPRC/SOEST, University of Hawaii) from the public site 
of NASA JPL PO DAAC and interpolated (linearly) onto exact locations and times of drifters. 
Optimal coefficient and angle of vector T to the QuikSCAT wind vector are shown in Figure 11, 
which reveals significant differences in both parameters from Figure 9. Taking into account that 
the QuikSCAT dataset starts in 1999, temporal distribution of the data is quite different from the 
one described in Section 2 that can partially be the source of discrepancies between the two 
statistics. The other reason for the difference can be the difference between the two wind 
products. To eliminate the effect of temporal bias, we compared QuikSCAT winds with the 
contemporaneous NCAR/NCEP wind vectors and results are shown in Figure 12. Correlation 
coefficient (Fig.l2a) mainly varies between 0.2 and 0.8 and higher for higher wind speeds, 
probably defined by the relative level of noise. R.m.s. deviation (Fig. 12b) between the two winds 
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is in the range from 2 to 8 m/s and exceeds the formal errors of the QuikSCAT winds, which are 
announced to be 2 m/s in the wind speed and 20® in the wind direction. Ratio of the r.m.s. 
deviation to the wind speed in the selected latitude-wind-speed bin is shown is Fig. 12c and 0.5 is 
its rather typical value. Angle (Fig.l2d) between the wind vectors is typically ± 2-5®, but at some 
latitudes and wind speeds exceeds ± 10®. Coefficient between the winds (Fig.l2e) is generally 
larger than 0.6 and less than 1 (QuikSCAT winds are stronger). Although Figures 12d and 12e 
indicate differences between the wind scale and direction comparable to the formal errors of the 
QuikSCAT dataset, they in fact reveal remarkable systematic differences between the two wind 
products. Further comprehensive study is required to understand how these systematic 
differences are distributed along the spatiotemporal spectra, what fraction of the error is due to 
the differences in physical models used to describe the atmospheric boundary layer and whether 
these errors can be reduced. 


4. Biases in the statistics based on combined datasets having different spatiotemporal 
resolution 

Accuracy of the momentum balance described in Section 1 is affected by the differences between 
spectral properties of three absolutely independent datasets (drifters, satellite altimetry and 
winds) used jointly. Correction of the imbalance with underthought method produces a danger of 
biases, some of which will be considered in this Section. 

Drifters are proved to provide measurements of the surface velocity at good accuracy. However, 
their Lagrangian data are distributed very irregularly in space and time. It is very difficult to get 
quasi-continuous drifter observations at specific Eulerian locations. In addition, oscillations on 
periods close to tidal and inertial are usually very energetic in the velocity data. Inertial 
oscillations are usually caused by sharp variations of local wind during passages of atmospheric 
eddies or fronts. Their dynamics is described by the balance between terms (I) and (II) (and in 
some cases (IV)) of the equation (2), while term (HI) is assumed to be small. In practice, inertial 
oscillations are not well-resolved by the satellite ARGO system determining drifter coordinates a 
few times every day (in some experiments, like WOCE, many drifter transmitters were on only 
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every third day). The kriging (one-dimentional optimal interpolation) used at NOAA/AOML to 
produce the 6-hourly dataset does not necessarily provide satisfactory description of the highr 
frequency velocity signal. Modem satellite altimeter supplies periodic data of the sea level 
anomaly at pretty amazing accuracy (few cm). However, satellite tracks are sparse in time and/or 
in space. Moreover, along-track data contain only one component of the vector (HI) in (2). To 
make the altimetry data more friendly and practical, the Aviso (1996) fulfilled objective mapping 
of merged multi-satellite observations onto regular spatiotemporal grid. The price paid for this 
convenience is too large minimal temporal (10-15 days) and spatial (100-300 km) scales 
resolved. 

Comparison between drifter and satellite altimetry data (Niiler et al., 2003a) shows that Aviso 
maps contain most of strong mesoscale eddies, but underestimate swirling velocities of the 
eddies. Niiler et al. (2003a) suggested to correct this difference by rescaling geostrophic velocity 
estimated from altimetry by regressing it linearly to the low-passed drifter velocities. Typically, 
the coefficient varied between 1.3 and 2 and the correction improved eddy statistics significantly. 
At the same time, large (185 cm) mean sea level difference across the Kuroshio Extension / 
Subarctic Front system, significantly exceeding expected value of 140-150 cm was recognized as 
an error of the method. Careful consideration did reveal the possibility of the bias occurring with 
the abovementioned correction. The problem with the correction procedure is that the same 
rescaling is applied to all spatial scales, while actual errors are on the smallest scale only. It is 
especially hard to suggest any reasonably simple improvement in the situation when velocity 
anomalies of different signs have different spatial scales. Figure 13 showing spatial distribution 
of the skewness of temporal variations of zonal component of geostrophic velocity proves that 
such a situation takes place in the reality. Multiple sources of the skewness include prevailing 
longitudinal variations of the velocity in the jets and preferred routs that eddies of different signs 
take in the ocean. However, it was estimated that errors due to such a bias are less than 10% and 
are not sufficient to explain the overestimations by Niiler et al. (2003a). 

In the course of further study, more serious and more complex bias has been recognized. The 
source of this sort of the bias is in the mixed Eulerian-Lagrangian statitistical methods used. 
Namely, as high frequencies are missing in the mapped Aviso altimetry, Niiler et al. (2003a) 
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used the averaging in spatiotemporal bins to replace 6-hourly along-drifter-track data by the local 
monthly mean “independent observations” in regular spatial bins. As the faster drifter is moving 
the more bins it passes during the same time period, relative number of the resultant 
“independent observations” with larger velocities is larger than the relative number of the 
original 6-hourly fixes with the same large velocities. It was also shown that additional necessary 
conditions for this bias are low density of drifters and underestimated geostrophic velocities. 
Detailed description of the this study is provided by Maximenko (2003), who showed that all the 
mentioned conditions are satisfied and the bias gives rise to 30-50% overestimation of the mean 
velocities of all main jet currents and sea level differences across them. 
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Figure 1. 1992-2002 absolute mean sea level (Niiler et al., 2003b). Contour interval is 10 cm. 
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Figure 3. Zonal mean -curl(3T / dz) /(Pog) as expected from (1-2) (gray strips) and 
estimated from the NCAR/NCEP wind and Ralph and Niiler (1999) formula (3) in 
the Pacific (a), Atlantic (b) and Indian (c) oceans. Units are 10'*^ s'\ 




GSFC00.1 MSS-GGM01S, cm (From Brian Beckley) 
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Figure 4. Mean sea level based on the first data of the GRACE mission computed at the NASA Goddard 
Snace Flipht Center to order and decree 90 Ccoiirtesv of Brian RecklevV Units are cm. 





Figure 5. Smoothed mean sea level computed at the NASA Jet Propulsion Laboratory (courtesy of Victor 
Zlotnicky). Also similar to filter mean sea level shown in Figure 3. Contour interval is 10 cm. 
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Figure 6. Difference between the sea levels shown in Figures 1 and 4. Contour interval is 10 cm. 
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Figure 7. Optimal angle and coefficient between dT / dz /(Po) and local NCAR/NCEP 
wind vector W as a function of latitude (dark blue). Red is for RN99 and green for drifters 

that lost their drogues. 
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Figure 9. Optimal angle and coefficient between 3X / dz /(Po) and local 
NCAR/NCEP wind vector W. Individual dots are for the 9® latitude x 3 m/s 
wind bins, centered on the 1° latitude x 1 m/s grid. Red are for the winds lower 
than 11 m/s. Black solid lines are RN99. 
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Figure 10. Correlation coefficient (CC) and optimal coefficient and angle 

between dT / dz /(Po) and local NCAR/NCEP 10m wind vector W. Patches 
represent 9° latitude x 3 m/s wind bins, centered on 1° latitude x 1 m/s grid. Only 

bins with CC > 0.1 are shown. 





























Skewness 


22 



.Skewness <U >/<U > of temporal variations of zonal component of geostrophic velocity estimated 

from the Aviso maps. 
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